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The  advanced  buildings  of  tomorrow  will  need  to  take  advantage  of  renewable,  ambient  and  waste 
energy  to  approach  ultra-low  energy  buildings.  Such  buildings  will  need  to  consider  Thermal  Energy 
Storage  (TES)  techniques  customized  for  smaller  loads. 

Recently,  TES  has  attracted  increasing  attention  due  to  the  potential  benefits  it  can  deliver  in  energy 
efficiency,  shift  load  from  peak  to  off-peak,  economics  and  environmental  impact.  Advanced  design 
tools  and  technical  improvements  are  required  in  TES  technology  and  systems.  Indeed  the  design  of  the 
building  and  the  TES  are  often  not  coordinated.  A  building  integrated  with  distributed  thermal  storage 
materials  could  shift  most  of  peak  load  to  off-peak  time  period. 

Even  though  various  tools  have  been  developed  to  model  the  behavior  of  applied  phase  change 
materials  (PCM),  unacceptable  accuracy  and/or  high  computational  time  are  addressed  as  their  major 
drawbacks.  This  implies  that  the  development  of  a  fast  and  reliable  model  is  necessary  in  order  to 
simulate  the  long-term  behavior  of  these  materials,  especially  for  design  and  optimization. 

Therefore,  a  new  and  fast  one-dimensional  analytical  model  is  proposed  in  this  paper.  The  PCM 
behavior  is  modeled  using  a  RC-circuit  concept  containing  variable  capacities  for  resistant  and 
capacitor.  In  this  approach,  the  length  of  mushy,  solid,  and  liquid  phases  in  each  period  of  time 
signifies  the  RC  capacity.  In  order  to  evaluate  the  performance  of  the  proposed  model,  a  computational 
fluid  dynamics  (CFD)  model  is  then  developed.  The  prediction  of  the  newly  developed  model  is 
compared  with  the  prediction  made  by  the  CFD.  There  are  good  agreements  between  the  predictions, 
and  the  results  clearly  show  the  high  performance  of  the  proposed  model. 
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1.  Introduction 

Building  sector  contributes  immensely  to  the  total  energy 
consumption,  particularly  for  its  space  conditioning  and  domestic 
hot  water.  In  Quebec,  70%  of  residential  buildings  use  electrical 


*  Corresponding  author.  Tel.:  +1  514  848  2424  93;  fax:  +1  514  848  7965. 
E-mail  address:  haghi@bcee.concordia.ca  (F.  Haghighat). 

1364-0321/$ -see  front  matter  ©  2012  Elsevier  Ltd.  All  rights  reserved. 
http://dx.doi.org/ 1 0.1 01 6/j.rser.201 2.04.053 


space  heating  system  which  accounts  for  30%  of  the  grid  critical 
peak,  while  91%  of  domestic  water  heaters  are  electric  and  they 
account  for  1750  MW  of  the  grid  peak  [1].  At  the  Canadian  level, 
the  domestic  water  heater  accounts  for  22%  of  the  total  energy 
use  for  households.  According  to  Hydro  Quebec  [2],  the  marginal 
cost  to  supply  electricity  during  the  peak  period  in  winter  is  10  $/ 
kW  and  will  be  increased  up  to  40  $/kW  by  2015.  At  the  same 
time,  in  Canada,  new  buildings  are  being  built  at  rate  of  1%  and 
old  buildings  stock  is  retrofitted  at  the  rate  of  2.2%  annually  which 
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Nomenclature 

L 

total  length  (m) 

H 

total  volumetric  enthalpy  (J) 

X 

space  coordinate  (m) 

Cp 

specific  heat  (J/kg  K) 

Subscripts 

t 

time  (s) 

T 

temperature  (K) 

S 

solid 

I< 

thermal  conductivity  (W/mK) 

F 

fluid 

a 

thermal  diffusivity  (m2/s) 

M 

mushy 

P 

liquid  fraction 

W 

wall 

2 

latent  heat  (kj/kg) 

P 

density  (kg  m-3) 

Non-dimensional  numbers 

Q 

heat  flux  (W/m2) 

E 

released  energy  (W/m2) 

Re 

Reynolds 

P 

dynamic  viscosity  (Pa  s) 

Nu 

Nusselt 

u 

Fluid  velocity  (m/s) 

Pr 

Prandtl 

s 

surface  position  (m) 

l 

length  in  one  specific  time  (m) 

means  by  2030  almost  50%  of  the  existing  would  be  retrofitted. 
Taking  this  information  into  account,  shifting  a  significant  portion 
or  the  whole  required  energy  for  space  heating  and  domestic 
water  heating  to  off-peak  periods  would  have  significant  eco¬ 
nomic  effects  on  both  energy  supply  and  demand.  This  shifting 
technique  is  accomplished  by  storing  energy  during  off-peak 
periods  to  be  utilized  during  peak  periods.  Building  envelope 
and  central  thermal  storage  have  been  used  as  thermal  energy 
storage  (TES).  Recent  study  also  showed  that  integration  of  PCM 
in  building  envelope  can  reduce  the  urban  heat  island  (UHI) 
phenomena  and  improve  pedestrian  comfort  and  health  [3,4]. 

Recently,  TES  has  attracted  increasing  attention  due  to  the 
potential  benefits  it  can  offer  in  energy  efficiency,  in  shifting  load 
from  peak  to  off-peak,  in  emergency  heating/cooling  load,  in 
economics  and  in  environmental  impact.  Nonetheless,  advanced 
design  tools  and  technical  improvements  are  required  in  TES 
technologies  and  systems.  Indeed  the  design  of  the  building  and 
TES  are  often  not  coordinated.  A  building  integrated  with  dis¬ 
tributed  thermal  storage  materials  could  shift  most  of  peak  load 
to  off-peak  time  period  as  it  is  important  to  plan  for  the 
requirements  of  the  buildings  of  the  future.  Hence,  research  is 
needed  to  ensure  that  the  actual  performance  of  buildings  with 
TES  in  use  matches  expectation  and  predictions  more  closely.  The 
use  of  TES  requires  smart  energy  management  in  buildings  to 
achieve  the  overall  goal  of  net  zero  energy  buildings.  It  is 
desirable  to  develop  technical  solutions  and  tools  that  incorporate 
advanced  TES  material  to  provide  alternative  option  to  conven¬ 
tional  solutions  and  to  meet  the  energy  demands  of  buildings  of 
the  future  ( SMART  NET-ZERO  BUILDINGS). 

As  part  of  Annex  23,  an  international  survey  was  carried  out  to 
review  the  actual  implementation  of  TES  use  in  energy  efficient 
buildings.  The  results  showed  that  the  sensible  thermal  storage  (e.g., 
water  heating  system  and  geothermal)  and  latent  thermal  storage 
(ice)  are  mature  technologies  with  mostly  existing  reliable  design 
and  system  integration  tools.  The  results  of  the  survey  also  indicated 
that  there  is  a  lack  of  knowledge  and  design  tools  for  sizing, 
optimization,  system  integration  and  control  of  other  type  of  latent 
thermal  storages  such  as  phase  change  materials  (PCM).  PCM  has 
been  considered  as  excellent  candidate  to  be  used  as  TES  due  to  its 
high  capacity  to  store  energy.  It  can  be  integrated  as  layer  in  a  wall 
[5,6],  in  a  floor  or  ceiling  combined  by  a  heating  system  [7]. 

Despite  many  publications,  only  very  limited  integration  of  PCM  in 
existing  buildings  can  be  found,  and  there  is  almost  no  information 
regarding  validated  design  tools  for  such  applications.  Most  research¬ 
ers  have  used  in-house  solution  without  any  systematic  comparison 
with  other  alternative  design  approaches  [8].  It  is  desirable  to  develop 


technical  solutions  and  tools  that  incorporate  advanced  TES  to 
provide  alternative  energy  option  to  conventional  solutions  and  to 
meet  the  energy  demands  of  future  buildings.  Also,  the  existing 
techniques  for  so  called  “passive  design”  require  to  be  pursued  with 
greater  vigor,  supported  by  improved  building  thermal  simulation  to 
optimize  their  performance.  PCM  can  be  embedded  within  walls  and 
other  construction  materials  to  moderate  temperature  variations  and 
act  to  enhance  thermal  capacity  as  well  as  to  shift  load  from  peak  to 
off-peak.  Application  of  PCM  in  building  construction  materials  has 
not  been  fully  exploited  and  is  generally  underutilized.  Further 
research  has  shown  that  the  application  of  PCM-impregnated  gypsum 
board  in  a  naturally  ventilated  passive  solar  building  reduced  the 
heating  demand  up  to  90%  during  autumn  and  spring  [9].  It  was 
concluded  that  further  research  is  needed  to  make  it  efficient  for 
heating  season. 

The  performance  of  the  passive  design  depends  on  the  func¬ 
tion  of  building,  office  or  residential  buildings,  and  its  relation 
between  the  PCM  type  and  its  position  with  the  building  envelope 
as  well  as  its  location  within  a  room  in  a  building.  For  example, 
PCM  can  be  either  introduced  in  the  insulation,  in  the  interior 
finishing  or  in  the  storage  (i.e.,  hot  water  tank);  further  research  is 
needed  to  optimize  the  whole  system.  Numerical  analysis  of 
thermal  effects  of  hollow  core  concrete  slabs  is  well  documented. 
Many  simulations  have  been  performed  to  determine  the  effects 
of  air  flow  rate,  hollow  core  diameter  and  inlet  air  temperature  on 
the  energy  saving  and  thermal  load  shift  properties  of  hollow  core 
slabs  [10].  In  addition,  PCMs  has  been  incorporated  in  the 
construction  of  hollow  core  slabs.  Moreover,  their  effects  on  the 
slab  energy  saving  in  a  mild  climate  have  been  investigated  [11]. 
The  combination  of  thermal  mass  and  phase  change  materials 
shows  great  promise  in  terms  of  energy  saving.  Efficient  simula¬ 
tion  tools  are  needed  to  investigate  the  range  of  parameters  and 
the  impact  different  climate  on  its  performance,  and  the  cost. 

Analysis  and  design  of  TES  for  future  buildings  has  to  account 
for  several  interactions  of  complex  physical  processes  inherent  to 
this  system.  Although  existing  building  simulation  tools  can  model 
building  envelopes,  thermal  mass  or  duct  systems,  they  cannot 
efficiently  and  accurately  model  double  skin  facades  integrated 
with  TES  [12]  and/or  dynamic  insulation  walls  imbedded  with 
PCM  [13].  These  tools  cannot  also  accurately  model  the  complex 
interaction  processes  between  PCM  and  the  building. 

Computational  Fluid  Dynamics  (CFD)  codes  have  been  widely 
used  in  the  simulation  of  PCM.  Experimental  data  have  been  used  for 
qualitative  and/or  quantitative  validation  of  the  CFD  simulations  [8]. 
Though  the  prediction  by  these  models  is  generally  similar  to  that  of 
CFD,  they  are  not  efficient  design  tools.  Development  of  a  validated 
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simplified  model  is  a  challenging  task  since  modeling  PCM  involves 
transient  non-linear  phenomena  with  a  moving  liquid-solid  interface 
whose  position  is  unknown  a  priori.  The  existing  approaches  are  not 
efficient  for  the  design  of  an  energy  efficient  building  integrated  with 
TES.  Most  of  existing  PCM  models  for  integration  in  buildings  and 
even  those  being  implemented  in  building  simulation  programs 
(TRNSYS,  EnergyPlus)  present  some  major  limitations.  As  an  example, 
the  model  developed  by  Ibanez  et  al.  requires  experimental  data  for 
the  system  to  be  designed,  which  cannot  be  readily  used  for  new  and 
innovative  designs  and  optimization  [14].  Also,  the  3D  model  devel¬ 
oped  by  Lamberg  et  al.  requires  extensive  calculation  time,  which 
makes  its  use  impractical  for  design  and  optimization  purposes 
[15,16].  Dutil  et.  al  carried  out  a  comprehensive  review  of  the  existing 
PCM  models  and  reported  that  “many  of  the  recent  studies  discuss  their 
results  qualitatively  only  as  the  comparison  with  a  graph  taken  from 
publications...”  [8].  Therefore,  a  new  generation  of  efficient  design 
tools  and  their  comprehensive  validation  are  needed. 

This  paper  reports  the  development  of  a  simplified  model  for 
application  of  PCM  TES  in  the  building  simulation  programs.  It  is  a 
fast,  efficient  and  validated  model  that  can  be  used  in  the  design 
of  an  energy  efficient  building  or  TES. 


2.  Brief  review  on  existing  models 


Stefan  problem  is  attributed  to  type  of  problems  with  chan¬ 
ging  phase  boundary  with  time.  Various  analytical  solutions  in 
different  situation  were  first  introduced  by  Jozef  Stefan  [8,17]. 
These  solutions  are  mostly  one-dimensional  and  defined  in  semi¬ 
infinite  or  infinite  regions.  The  assigned  boundary  and  initial 
conditions  are  also  simple  and  mostly  constant  while  the  mushy 
zone  is  not  considered  in  solutions. 

Numerical  techniques  have  been  carried  out  to  model  not  only 
the  moving  boundary  of  phases,  but  also  to  calculate  the  velocity 
of  mushy  region.  In  one  approach,  the  main  assumption  is  Stefan 
condition  in  which  latent  heat  either  can  be  absorbed  or  desorbed 
in  boundaries: 


STs 

St 


I< 


ST 


f 


f 


St 


However,  it  should  be  noted  that  discontinuity  could  occur  in 
numerical  model  between  solid  and  fluid  phases  due  to  a  differ¬ 
ence  in  physical  properties.  As  another  technique,  the  behavior  of 
mushy  zone  can  be  simulated  using  enthalpy  approach  [18].  In 
this  model,  conservation  of  energy  can  be  expressed  in  terms  of 
total  volumetric  enthalpy  and  temperature  as  bellow: 

H=  [  CpdT+pfX  (2) 

Jrre{ 


Here,  the  first  and  second  parts  represent  the  sensible 
enthalpy  and  the  latent  heat,  respectively.  Therefore,  the  energy 
equation  in  the  PCM  can  be  written  as: 

^(pH)  +  V.(p~uH)  =  V.(KVT)+S  (3) 

where  S  is  the  source  term  and  different  treatments  have  been 
proposed  to  address  this  term  specially  in  mushy  region. 

Different  techniques  (i.e.,  finite  difference,  finite  element,  and 
finite  volume)  are  also  widely  used  to  solve  enthalpy-based 
equations  of  PCM.  Even  though  the  results  are  mostly  acceptable 
and  reliable  for  a  wide  range  of  engineering  problems,  dealing 
with  above  mentioned  methods  is  not  always  a  straight  forward 
procedure,  and  oscillation  of  temperature  in  some  occasion  is 
broadly  reported  [19].  Nonetheless,  advanced  computational 
techniques  (e.g.,  using  moving  meshes  or  refining  meshes  within 
boundary  regions)  helped  to  improve  the  accuracy  of  these 
methods.  The  main  drawback  of  such  techniques  is  still  related 


to  their  intensive  time  and  CPU  cost,  where  annual  simulation  is 
mainly  deemed  for  building  energy  calculation. 

Furthermore,  other  forms  of  the  continuity  and  momentum 
equations  are  adapted  to  model  the  PCM.  For  example,  the 
temperature  transforming  model  (TTM)  keeps  velocity  at  zero 
[20]  or  variable  viscosity  of  the  medium  (VVM)  approach  sets  the 
viscosity  within  a  high  value  in  the  solid  phase  [21].  Both 
methods,  however,  again  have  difficulties  in  conserving  continu¬ 
ity  and  code  development. 

Transfer  function  is  a  traditional  way  to  model  transient  effect  in 
building  energy  load  calculation  [22,23].  Recently,  a  conduction 
transfer  function  (CTF)  model  is  adapted  for  application  of  the  PCM 
in  buildings  using  response  factor  to  relate  the  current  and  past  heat 
flux  to  current  and  past  temperatures  of  the  building  element.  This 
procedure  assumes  that  materials’  density,  thermal  conductivity  and 
specific  heat  are  constant  for  the  period  of  a  simulation  where  this  is 
not  a  valid  assumption  for  the  PCM,  since  their  phase  and  properties 
are  continuously  changing.  Basic  formulation  of  this  method  pre¬ 
sumes  a  simple  relation  between  heat  capacity  and  temperature. 
The  accuracy  of  such  models  is  not  acceptable  although  it  might  be 
improved  assuming  various  layers  for  the  PCM.  Again,  shortcoming 
of  this  method  is  addressed  for  annual  thermal  load  calculations  and 
uncertainty  on  accuracy  improvement.  Energy  Plus,  DOE-2  and 
BLAST  are  comprehensive  and  computationally  efficient  tools  for 
annual  building  energy  load  calculations  which  are  developed  based 
on  the  CTF  method.  Barbour  and  Hittle  [24]  suggested  obtaining  a 
set  of  CTFs  and  switching  them  between  the  phases.  Their  results, 
however,  show  an  error  of  20%,  which  is  not  acceptable  for  the 
design  of  energy  efficient  buildings. 

Another  concept  is  to  model  the  PCM  with  one  or  various  set  of 
constant  resistant  and  capacitors  (RC)  even  though  assuming  one  RC 
circuit  cannot  provide  proper  accuracy.  On  the  other  hand,  more 
advanced  techniques  should  be  integrated  to  the  multiple  sets  of 
RCs’  model  (e.g.,  genetic  algorithm)  in  order  to  optimize  series  of 
capacitance  and  resistance  [25].  This  indicates  that  intensive  efforts 
are  needed  for  the  development  of  such  model  for  building  energy 
simulation  purposes.  Moreover,  expecting  higher  accuracies  requires 
increasing  the  number  of  RCs,  which  enormously  affect  the  compu¬ 
tational  efforts.  Table  1  briefly  summarizes  all  existing  models  for 
the  PCM  modeling  in  term  of  accuracy  and  computational  time  in 
building  applications.  Unlike  presented  models  in  Table  1,  the 
proposed  model  in  the  next  section  claims  to  provide  an  acceptable 
accuracy  and  a  low  computational  load. 


3.  Proposed  model 

As  discussed  earlier,  RC  models  are  reliable  techniques  even 
though  their  accuracy  significantly  depends  on  the  number  of 
synchronized  RC-circuits.  This  means  that  higher  number  of  RC- 
circuits  increases  simultaneously  their  accuracy  and  CPU  time. 
Therefore,  having  only  one  RC-circuit  would  not  satisfy  the 
expected  accuracy.  In  this  study,  a  new  one-circuit  RC  model  is 
developed  with  a  considerable  higher  accuracy.  Unlike  the  tradi¬ 
tional  models,  the  capacity  of  the  capacitor  and  resistant  varies 
with  time  in  the  current  model.  Depending  on  existence/ 


Table  1 

Available  approaches  to  model  phase  change  materials  in  building  applications. 


Approach 

Accuracy 

Time  and  CPU  cost 

Enthalpy  models 

Good 

Very  high 

TTM  and  WM 

Acceptable 

Very  high 

RC  models 

Fair 

High 

CTF 

Fair 

Medium 

Stefan  models 

Poor 

Low 
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Fig.  1.  Reduction  of  latent  heat  energy  during  solidification  process,  (a)  Left  wall:  S2  and  right  wall  S2mechanisms  and  (b)  Left  wall:  SI  &  S2  and  right  wall  S2  &  S3mechanisms. 


coexistence  of  mushy,  solid,  or  liquid  phases  in  each  time  step,  a 
value  can  be  assigned  to  a  single  RC-circuit.  For  example,  when 
solid  phase  changes  to  mushy  and  liquid  phases,  the  capacity  of  the 
RC-circuit  alters  accordingly.  The  length  of  each  phase  therefore 
presents  the  capacity  of  the  RC-circuit  during  each  period  of  time. 


3.1.  Solidification  process 


The  proposed  model  consists  of  three  phases;  fluid,  mushy, 
and  solid.  The  temperature  of  mushy  zone  is  assumed  to  vary 
linearly  between  solidification  (Ts)  and  melting  temperatures  (7}) 
of  the  PCM.  Therefore,  the  liquid  fraction  can  be  defined  as  below: 

P  =  (T-Ts)/(Tf-Ts)  (4) 


where  in  solid  phase  p=0  and  in  fluid  phase  ft  =  1 

While  it  is  assumed  that  thermal  stratification  can  be  neglected 
due  to  the  thickness  of  the  PCM,  the  heat  is  expected  to  one- 
dimensionally  conduct  through  the  wall  and  the  PCM.  As  depicted  in 
Fig.  1,  the  area  under  diagram  ft-x  shows  the  potential  energy 
deposited  in  mushy-phase  or  liquid-phase  of  the  PCM,  and  the 
alteration  of  this  area  thus  denotes  the  released  energy: 


A  E  = 


pXftdx 


Illustrated  in  Fig.  1(a)  and  (b),  the  PCM  with  a  greater  tempera¬ 
ture  than  Tf  has  initially  a  latent  energy  E=IpCMxpx/l(^l) 
During  solidification  process,  the  PCM  loses  part  of  its  energy  until 
it  reaches  its  solidification  temperature  (/?=1).  It  should  be  noted 
that  losing  energy  could  occur  either  by  a  decrease  in  level  of  liquid 
fraction  or  an  increase  in  length  of  mushy-zone.  In  general,  the 
conducted  heat  fluxes  through  the  left  and  right  walls  define  the 
form  of  the  energy  release/gain,  which  can  be  in  a  combination  of 
the  following  three  mechanisms: 


51 - Decrease  in  level  of  the  liquid  fraction  in  mushy-phase 

52- Conversion  of  liquid-phase  to  mushy-phase 

53- Conversion  of  mushy-phase  to  solid-phase 

Since  liquid  fraction  (/?)  changes  linearly  between  solid-phase 
and  liquid-phase,  the  energy  release  defined  in  Eq.  (5)  can  be 
shown  by  different  triangles  in  Fig.  1.  For  example,  triangles  1  and 

2  (Fig.  la)  represent  released  energy  due  to  change  of  liquid- 
phase  to  mushy-phase  (S2),  and  as  depicted  in  Fig.  1(b),  triangles 

3  and  5  are  again  related  to  alteration  of  the  liquid-phase  to 
mushy-phase.  On  the  other  hand,  triangle  4  shows  how  energy  is 
released  when  level  of  the  liquid  fraction  is  decreased  in  the 
mushy-phase  (SI).  Moreover,  triangle  6  displays  released  energy 
related  to  change  of  the  mushy-phase  to  solid-phase  (S3). 

The  energy  release  process  similar  to  triangles  5  and  6  (right  wall 
in  Fig.  1(b))  is  due  to  phase  change  of  mushy  to  solid  and  liquid  to 
mushy.  As  depicted  in  Fig.  2,  therefore,  Emushy^soiid  can  be  calculated 
as  the  area  between  line  1  and  2  or  hatched  triangle.  Similarly, 
Eiiquid  — >  mushy  is  related  to  the  area  between  line  2  and  3  or  shaded 
triangle  (Fig.  2).  On  the  other  hand,  the  length  of  mushy  zone  is 
initially  Lm  =  /m  +  /s.  This  means  that  the  initial  profile  of  the  mushy 
region  can  be  exhibited  by  line  1.  Then,  part  of  the  mushy-PCM 
energy  (equal  to  hatched  triangle)  is  released.  Therefore,  line 
2  depicts  a  new  linear  profile  of  /?.  The  new  length  of  mushy  zone 
is  also  equal  to  lm.  Eventually,  another  portion  of  liquid-PCM  energy 
(equal  to  shaded  triangle)  is  released  which  changes  the  length  of 
the  mushy-PCM  to  lm+lf.  Line  3  illustrates  the  final  linear  profile  of  ft. 
Therefore,  the  whole  process  implies  that  the  level  of  mushy-PCM 
energy  is  reached  from  line  1  to  line  3.  The  same  analogy  can  be 
used  to  explain  the  behavior  of  a  process  similar  to  Fig.  1(b)  on  the 
left  wall  where  part  of  the  energy  release  is  due  to  the  change  of 
liquid-phase  to  mushy-phase  (triangle  3)  and  part  is  due  to  a 
decrease  in  the  level  of  the  liquid  fraction  (triangle  4). 
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3.2.  Melting  process 

The  same  analogy  can  be  used  to  model  melting  process  of  the 
PCM.  Conducted  heat  transfer  through  walls  is  accumulated  in  the 
PCM  in  three  different  forms: 

Ml -Increase  in  level  of  the  liquid  fraction  in  mushy-phase 

M2-Conversion  of  mushy-phase  to  liquid-phase 

M3-Conversion  of  solid-phase  to  mushy-phase 

Obviously,  the  formation  of  the  above  mentioned  mechanisms 
or  their  combinations  depends  on  heat  fluxes  in  both  left  and 
right  walls.  For  example,  Fig.  3  demonstrates  a  situation  in  which, 
first,  both  walls  absorb  energy  through  mechanism  M3.  This 


Fig.  2.  Energy  release  with  length  change  of  the  mushy-PCM. 


implies  that  triangles  1  and  2  show  the  amount  of  energy 
accumulated  in  the  PCM  through  mechanism  M3.  Second,  Ml 
and  M3  mechanisms  contribute,  respectively  to  deposit  energy  in 
the  left-side  of  the  PCM  where  M2  and  M3  are  attributed  to 
absorb  energy  in  the  right-side  (Fig.  3(b)).  Again,  triangles  3  and 
6  represent  the  growth  of  the  mushy  region.  Triangle  4,  however, 
is  related  to  Ml  or  the  increase  in  the  level  of  liquid  fraction  in 
mushy-phase.  Eventually,  triangle  5  depicts  the  absorbed  energy 
due  to  transformation  of  the  solid-phase  to  mushy-phase. 

3.3.  Methodology 

Modeling  of  the  PCM  behavior  with  variable  capacitor  concept 
is  explained  in  this  section.  First,  in  the  proposed  model,  each 
phase  is  assigned  to  one  node.  For  example,  Fig.  1(a)  and  (b)  can 
be  simulated  with  two  and  three  nodes,  respectively,  or 
Fig.  3(a)  and  (b)  can  be  presented  with  two  and  three  nodes  since 
solid-mushy  and  solid-mushy-liquid  phases  coexist  together, 
respectively.  Then,  energy  balance  can  be  written  for  each  node 
assuming  variable  length  for  existing  phases.  The  variable  lengths 
are  later  considered  to  calculate  the  capacity  of  the  RC-circuit  in 
each  period  of  time. 

For  example,  using  electrical  network  method,  the  proposed 
one-dimensional  model  for  solidification  process,  when  all  three 
phases  coexist  together,  is  illustrated  in  Fig.  4.  Presented  in 
Eq.  (6),  a  portion  (Zs)  of  solidified  mushy-PCM  (q2)  in  addition  to 
conducted  heat  from  mushy-zone  (g3)  supply  conducted  heat 
through  the  left  boundary  condition  (g!  =  qib).  Exhibited  in  Eq.  (7), 
the  conducted  heat  from  mushy-zone  (g3)  is  simultaneously 


Fig.  3.  Accumulation  of  latent  heat  energy  during  melting  process,  (a)  Left  wall:  M3  and  right  wall  M3  mechanisms  @  time  t  and  (b)  Left  wall:  Ml  &  M3  and  right  wall  M2 
&  M3  mechanisms  @  time  t  +  dt. 
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provided  by  alteration  of  a  portion  (If)  of  liquid-PCM  to  mushy- 
PCM  (q4),  and  transfer  of  generated  heat  (q6),  if  the  liquid-PCM 
has  higher  temperature  than  7},  through  the  mushy  zone  inte¬ 
grated  with  coming  heat  from  right  boundary  condition  ( qrb ). 
Eventually,  conversion  of  phases’  length  to  each  other  has  to  be 
conserved  according  to  Eq.  (9). 


Q2  +  ^3  —  —  Qlb  => 


( Tf-Ts ) 

Lm  +  If— Is 
Km 


+  E 


mushy  ->  solid 


(Ts-Tw) 

Ls  +  Is 

Ks 


Q5+Q4  ~  Q3  => 


T-Tf 

k 

Nu.I<f 


liquid. 


>  mushy 


(Tf-Ts) 

Lm  +  If —Is 
Km 


P3  ~  Ps  + 


(Pf-Ps) 
lm  +  If 


(X~ls), 


X  Is 

lm  +  If 


(13) 


Hence,  the  alteration  of  energy  in  mushy-PCM  through  the 
hatched  and  shaded  triangles  can  be  obtained  from  the  following 
equations: 


2 


Q4  —  E  liquid  -» mushy  —  ^ 


‘h  +  In 


(P2P2—  P3ft)d*  + 


’  h  +  Li  +  If 
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( Pf-PsP3 )dx 

(14) 


Q3  —  Emushy^>  solid  — 
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(PlPl-P2p2)dx 


(15) 


Qe+Qrb  —  <?5  =*“  PfCpfLfT  +  c7r6  — 


T-Tf 

k 

Nu.I(f 


Epcm  =  Lm  J-Ls+Lf  (9) 

where  LPC m  is  the  total  length  of  the  PCM  and  is  a  constant 
number.  The  following  empirical  equation  can  be  also  used  to 
express  Nusselt  number  [26]: 


Again,  Eqs.  (14)  and  (15)  are  developed  based  on  unknown 
variables  of  Zm,  ls ,  and  lf.  Thus,  Eqs.  (5)  through  (15)  can  be  used  to 
find  these  variables  at  each  time  step  with  calculation  of  new  IM, 
Ls  and  Lf  as  below: 


(Es) new  —  (Es)old  Is 

(16) 

(Ef  )new  =  (Ef  )o\d~  If 

(17) 

Nu  =  0.664  x  Re0  5  x  Pr°  33 


(10) 


As  shown  in  Eqs.  (6)  through  (10),  left  and  right  boundary 
conditions  can  be  easily  assigned  to  this  model  in  order  to  find  the 
unknown  variables;  ls,  lf,  and  T.  However,  to  solve  these  equations, 
it  is  necessary  to  first  obtain  the  released  energy  from  phase 
changes  within  the  PCM,  including  Emushy^solid  and  Eiiquid^mushy. 

To  calculate  AE  presented  in  Eq.  (5),  one  can  assume  linear 
regression  for  changing  density  and  liquid  fraction  of  the  PCM  in 
mushy  zone  as  bellow: 


P 1  —  Ps  + 


(Pf-Ps) 
Is  +  lm 


X 

Is  +  lm 


(11) 


,  (Pf 

P2  ~  Ps  ^  j 


■Ps) 


(x-lsl  p2  = 


X-L 


lm 


(12) 


Fig.  4.  Schematic  of  proposed  1 -dimensional  model. 


(Em) new  —  (Em) old  Is  If 


(18) 


4.  Results 

The  performance  of  the  proposed  one-dimensional  model  is 
validated  with  a  benchmark  (task  C)  defined  by  ANNNEX  23  [27]. 
This  task  includes  several  one-dimensional  case  studies  in  order  to 
simulate  the  behavior  of  a  wallboard  containing  the  PCM.  In  test  cases 
of  the  task  C,  the  long-wave  and  short-wave  radiations  are  assumed 
to  be  neglected.  Moreover,  the  convective  heat  transfer  coefficients 
are  assigned  as  two  constant  values  both  for  left-side  and  right-side 
boundary  conditions  (Table  2).  Here,  the  case  study  has  selected  a 
10  mm  sheet  of  the  PCM  without  insulation  (case  study  P10).  PCM 
characteristics  in  addition  to  initial  and  boundary  condition  are  also 
provided  in  Table  2  (case  A).  It  should  be  mentioned  that  the 
proposed  model  is  indirectly  validated  with  P10  case  study  as 
demonstrated  in  Fig.  5.  It  means  that  a  two-dimensional  computa¬ 
tional  fluid  dynamics  (CFD)  model  based  on  enthalpy  equations  is 
first  implemented  to  model  the  above  mentioned  case  study.  The 
main  reason  for  carrying  out  the  indirect  validation  is  for  comparing 
the  speed  of  the  proposed  model  with  a  two-dimensional  accurate 
but  time  consuming  model. 

A  structured  40  x  500  mesh  domain  (0.01  m  x  0.21  m)  is  gen¬ 
erated  for  two-dimensional  CFD  domain.  The  buoyancy  effect  is 
not  considered  in  order  to  convert  the  problem  to  a  one-dimen¬ 
sional  case.  Moreover,  the  conductivity  is  presented  with  two 


Table  2 

PCM  characteristics,  initial,  and  boundary  conditions. 


The  PCM  characteristic  and  Boundary  Condition 

TASK  C 

Case  A 

Solidification 

Case  B 

Melting 
Case  C 

Length  (m) 

0.01 

0.01 

0.01 

Left-side  convective  heat  transfer  coefficient  (W/m2I<) 

2.5 

10 

5 

Right-side  convective  heat  transfer  coefficient  (W/m2I<) 

8.0 

10 

20 

Density  (kg.m2) 

1100 

1100 

860 

Specific  heat  (J/kg  K) 

80,000 

140,000 

100,000 

Conductivity  (W/mK) 

0.25  (T<  285) 

0.25  (T<  285) 

0.23 

0.20  (285  <  T) 

0.20  (285  <  T) 

0.23 

Solidification  temperature  (K) 

281 

287 

290 

Melting  temperature  (K) 

301 

289 

299 

Left-side  air  temperature  (K) 

293 

283 

300 

Right-side  air  temperature  (K) 

285  +  20t/3600 

283 

302 

The  initial  PCM  temperature  (K) 

285 

290 

285 
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. . . . 


C 


. 


Fig.  5.  One  dimensional  modeling  of  the  PCM  sheet  (P10)  with  coexistence  of 
(a)  only  one  phase,  (b)  two-phases  and  (c)  three-phases. 


numbers  below  and  above  285  K.  Second-order  upwind  is  also 
considered  as  discretization  scheme  for  the  momentum.  Further¬ 
more,  a  convergence  criterion  is  applied  to  keep  residuals  less 
than  10“ 5  for  energy  and  continuity  equations. 

A  MATLAB  code  is  developed  to  evaluate  the  performance  of 
the  proposed  model.  As  discussed  earlier,  the  number  of  nodes 
represents  the  coexistence  of  phases  (Fig.  5).  When  the  PCM  only 
contains  one  phase  (pure  solid  or  liquid),  one  node  is  used  to 
simulate  its  behavior.  With  a  same  analogy,  two  or  three  nodes 
are  employed  when  two  (solid-mushy  or  liquid-mushy)  or  three 
phases  (solid-mushy-liquid)  coexist  together.  The  distance 
between  nodes  is  also  calculated  based  on  equation  given  in  the 
previous  section.  The  proposed  model  automatically  switches 
between  one,  two,  and  three  nodes  depending  on  the  temperature 
imposed  by  left  and  right  boundaries.  In  other  words,  the  length 
of  each  phase  provides  the  existence  of  one,  two,  or  three  nodes 
and  the  required  distance  between  them. 

The  validated  CFD  model  is  later  used  for  evaluation  of  the 
proposed  model.  For  this  purpose,  two  more  cases  (B  and  C)  are, 
respectively  defined  for  solidification  and  melting  to  show  the 
performance  and  flexibility  of  the  proposed  model. 

As  illustrated  in  Fig.  6,  the  results  of  CFD  model  for  left-wall 
temperature  show  good  agreement  with  corresponding  tempera¬ 
ture  of  the  above  mentioned  case  study  of  task  C.  The  average 
discrepancy  of  the  CFD  model  from  provided  results  is  less  than 
0.4  K  with  maximum  error  of  0.8  K.  After  validation  of  the  CFD 
model  with  the  outcomes  of  Task  C,  the  CFD  model  itself  was 
applied  for  evaluation  of  the  proposed  model. 

Solidification  process  is  simulated  with  both  CFD  and  proposed 
models  (Case  B).  As  shown  in  Table  2,  the  case  study  B  is  similarly 
defined  for  one  layer  PCM  wall.  The  characteristic  of  the  PCM  and 
boundary  conditions  are  completely  different  from  the  validation 
case  study  of  task  C  in  order  to  exhibit  the  flexibility  of  the 
proposed  model. 

The  comparison  of  two  models  is  depicted  in  Fig.  7.  It  can  be  seen 
that  the  proposed  model  perfectly  simulates  the  temperatures  above 
solidification  temperature  (TS=297I<).  The  mean  and  maximum 
errors  are,  respectively  around  0.11  K  and  0.60  K  during  solidifica¬ 
tion  process.  When  temperature  fluctuates  between  7}  and  Ts,  the 
Case  B  has  discrepancy  below  one  percent  except  for  small  numbers 
close  to  the  solidification  temperature  (Ts).  This  error  is  related  to 
the  change  of  entire  mushy-phase  to  solid-phase,  which  causes  a 
sudden  decrease  in  resistance  through  the  PCM.  The  trend  of  the 
temperatures  below  Ts  and  duration  of  whole  process  are  also 
successfully  simulated  with  the  proposed  model  (Fig.  7). 

Furthermore,  the  performance  of  the  proposed  model  to 
simulate  the  melting  process  is  examined  by  introducing  a  new 
case  study  (C)  with  a  new  PCM  as  shown  in  Table  2.  Here,  the 
solidification  and  melting  temperatures  are  selected  to  be  290  I< 


Time  (s) 

Fig.  7.  Comparison  of  the  proposed  (Case  B)  and  CFD  models  for  solidification 
process. 

and  299  K,  respectively.  As  exhibited  in  Fig.  8,  both  CFD  and 
developed  models  (Case  C)  present  the  melting  of  the  PCM  in 
three  parts:  below  Ts,  between  Ts  and  7),  and  above  7). 

The  proposed  model  can  fairly  follow  the  obtained  tempera¬ 
tures  in  left-wall  (Fig.  5)  of  the  CFD  model  since  the  turning  point 
in  three  main  parts  are  almost  the  same  for  both  models.  Both 
models  reach  Ts=297  I<  at  the  same  time,  although  the  Case  C 
reaches  7)  around  only  400  s  prior  to  the  CFD  model.  The  greatest 
discrepancy  is  related  to  the  point  where  the  temperature 
elevates  from  7}=  299  K.  The  mean  and  maximum  errors  are  also 
calculated  as  0.44  I<  and  1.30  K,  respectively. 


5.  Conclusion 

The  thermal  energy  storage  must  be  designed  to  respond  to 
the  actual  need  of  the  building.  To  design  a  TES  that  could  meet 
the  building  of  the  future  requirements,  it  must  be  an  integral 
part  of  the  system  of  the  building  design.  This  requires  a  reliable 
and  fast  simulation  tool  that  can  be  used  for  design  and  optimiza¬ 
tion  application. 

A  new  and  fast  one-dimensional  analytical  model  is  proposed 
for  simulating  PCM  behavior  using  a  RC-circuit  concept  contain¬ 
ing  variable  capacities  for  resistant  and  capacitor.  In  this  method, 
the  length  of  mushy,  solid,  and  liquid  phases  in  each  period  of 
time  signifies  the  RC  capacity.  In  order  to  evaluate  the  perfor¬ 
mance  of  the  proposed  model,  a  CFD  model  is  first  developed  and 
validated  with  a  well-defined  benchmark  [27];  task  (C).  After 
successful  development  of  the  CFD  model,  the  solidification  and 
melting  process  of  the  proposed  model  are  validated  under 
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various  PCM  characteristics,  boundary,  and  initial  conditions.  The 
obtained  results  show  the  considerable  capability  of  the  proposed 
model  in  order  to  simulate  the  behavior  of  the  PCM.  The  obtained 
mean  errors  for  solidification  and  melting  process  are,  respec¬ 
tively  0.11  K  and  0.44  K. 

There  are  good  agreements  between  the  predictions,  and  the 
results  clearly  show  the  high  performance  of  the  proposed  model. 
Unlike  simple  RC-circuit  models,  the  advantage  of  this  model  is 
due  to  its  variable  resistance  and  capacity  assigned  as  elements  of 
the  PCM.  Furthermore,  expanding  the  proposed  ID  model  to  a  2D 
space  would  be  a  ground  breaking  research  topic. 
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